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Abstract 



A new algorithm is presented for computing a direct solution to a system of 
consistent linear equations. It produces a minimum norm particular solution, 
a generalized inverse (of type {124}), and a null space projection operator. In 
addition, the algorithm permits an online formulation so that computations 
may proceed as the data become available. The algorithm does not require 
the solution of a triangular system of equations, nor does it rely on block 

^j. ■ partitioned matrices. 
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1 Introduction 

The problem of solving a system of linear equations is widespread across 
mathematics, science and engineering. When there are m equations in n 
unknowns, the equations are written as Ax = b, where A = (a#) is the 
coefficient matrix, x = (xj) is the vector of unknowns, and b = (bi). (The 
indices are i = 1, . . . , m and j = 1, . . . , n.) The entries in A, x and b are 
in a field J 7 , which is assumed to be endowed with an inner product. The 
examples in this paper assume J 7 = C. 
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2 THE ALGORITHM 



It is known that a solution to Ax = b may be written asz = x p +Xh, where 
the particular solution x p satisfies Ax p = b and the homogeneous solution Xh 
satisfies Axh = 0. Note that x^ is an element of the null space of A (i.e., 
Af(A)). When the dimension of M{A) is greater than zero or when A is 
rectangular, the inverse to A does not exist, and hence it's no longer possible 
to write the solution in the cogent form x = A~ 1 b. Nevertheless, when A~ x 
doesn't exist, one may still use the generalized inverse (G), which allows the 
two parts of the solution to be written as: 

x p = Gb 

x h = Py 
P = I n -GA 

where /„ is an n-by-n identity matrix, P is a null space projection operator 
and y is an arbitrary n-dimensional vector. (By definition, Py G M{A) 
for any y G J 7 ™.) Different properties of the solution Gb can be inferred 
by the type of the generalized inverse. (This is discussed in the literature 
[HI IU [21 Q] and will not be reviewed here.) The G that will be computed 
here is of type {124}, which among other things yields a solution Gb which 
has minimum Euclidean norm. 

There are other interesting features, such as the fact that the algorithm 
may be formulated in an online mode, which is defined to mean that the 
solution x can be formed while the input data (A and b) are still being 
acquired. This is useful in cases when it takes a relatively long time to 
access/create the input data. 

2 The algorithm 

The solution will require A to have its rows orthonormalized, perhaps using 
the usual Gram-Schmidt orthogonalization procedure. However, before dis- 
cussing a modified orthonormalization procedure, it is helpful to define the 
following. 

Definition: A quasi-orthonormal list of vectors consists of vectors with 
norm equal to 1 or 0. Those vectors of norm equal to 1 are mutually orthog- 
onal. 
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Definition: The row orthonormalization procedure (ROP) of a matrix 
Q is performed by left multplying it by a series of matrices M s (s = 1,2,3, ...) 
which affect an orthonormalization procedure on the nonzero rows of Q. The 
result is that the rows of Q form a quasi-orthonormal list of vectors. 

Thus, after applying the ROP, the rows are either zero or part of an or- 
thonormal set. (Pseudocode for the ROP may be found in Appendix A.) 
The ROP will here be applied to A as it exists in the context of the equa- 
tion Ax = b, causing the same operations to be applied to b as well. It is 
implemented with the nonsingular, m-by-m matrices M s which results in 

(■■■M 2 M 1 )Ax=(---M 2 M 1 )b (1) 

The Mi are applied to A, causing it to become A' = MA, where M = 
(• • • M 2 Mi). (Note that it is not required that the M s represent elementary 
row operations.) There are two versions for how to view the right-hand side 
of the previous equation: (1) apply the M s to b so that it becomes b' (where 
b' equals Mb); (2) accumulate the M s as a prefactor so that it becomes Mb. 
To summarize, after applying the ROP to Ax = b, these two variations may 
be written as 

A'x = b' (2) 

or 

A'x = Mb (3) 

The significance of the difference has mainly to do with computer implemen- 
tation issues, with respect to storage and reusability. 

Definition: Let Wq be a set which consists of the indices of the non-zero 
rows of an m-by-m matrix Q. Note that Wq C {1,2, ...,ra}. Define an 
index matrix by the following: 

(T).. = i 1 if i = i and i e Wqi 

[0 otherwise. 

This matrix is essentially just an identity matrix, except now the i-th row 
has a 1 only if i e Wq. This matrix is immediately applicable to A', whose 
rows comprise a quasi-orthonormal list of vectors. The following identities 
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may be easily verified: 

A'(A'y=l A , (4) 

A' = T A , A' (5) 

(Ay = (Ay z A , (6) 

Ia> V = b' (7) 

where the superscript * represents a conjugate transpose. 

Lemma: An arbitrary vector x G J 7 ™ may be expressed as 

x = (A')*w + x h (8) 

where A is an m-by-n matrix, w G J 7171 and Xh satisfies Axh = 0. 

Proof. Beginning with A'x = b', observe that A' : T n — > J 7171 is a linear 
transformation. It is then known (Thm. 18.3 [5]) 

T n = range[(A')*] © nul\[A'] (9) 

This direct sum decomposition may be used to rewrite an x G T n as 

x = w\ + w 2 (10) 

where 

Wi G range[(A')*] (11) 

w 2 G null [A'] (12) 

However, every vector that is in the range of some matrix Q may be expressed 
as Qv, for an appropriate v. Thus, for some w G J 7 " 1 , it follows that W\ may 
be written as 

Wl = (A')*w (13) 

Also, because M is nonsingular, nullLA'] = nullfA], so that w 2 G nullLA]. 
Substituting f lT3|) into f fTUj) gives the desired result. □ 

What remains is to determine w and to find a means of computing Xh- 
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Theorem: A solution to Ax = b is 

x = (A')*b' + x h (14) 

where all variables are as defined earlier. 

Proof. The proposed solution is verified by substitution into Ax = b 

Ax = A{A')*b' + Ax h 
= M^A'iA'Yb' 
= M' x T A ,b' 
= M- l b' 
= b 

D 

One could also take a constructive approach toward obtaining the above 
solution. Upon substituting (jSJ) into A'x = b' one obtains 

A'{A')*w = b'. 

which upon using (|3| becomes %&> w = b'. Equation ([Sj) may now be re- 
expressed as 

x = (A')* w + Xh 
= (A')*Z A > w + x h 
= (A>)*b' + x h 



Discussion 

One of the first things to notice is that (A')*b' is a particular solution to 
Ax = b, as well as a minimum norm solution. The reason for the latter is 
because it is an element of range[(yl / )*] which is the orthogonal complement 
of the null space of A' . In other words, the particular solution is orthogonal 
to the null space. The minimum norm nature of the particular solution will 
be revisited when the Penrose identities are checked. 

The first variation (Eqn. (J2J)) transforms Ax = b into A'x = V . In terms 
of augmented matrices, one applies the ROP to [A\b] to produce [-4'16'j, where 
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the rows of A' form a quasi-orthonormal list of vectors. The solution is formed 

as 

x v = (A')*bf 

Xh = Py 
P = I n - {Ay A' 

where y is an arbitrary vector in J 7 ™. This approach may be preferred in 
digital computations, when creating storage for M may be an issue. However, 
this variation would not be preferred if a solution is sought for additional b 
vectors, since the ROP step would have to be repeated. 

The second variation (Eqn. (J3J)) transforms Ax = b into A'x = Mb. 
In terms of augmented matrices, one applies the ROP to the augmented 
matrix [A|/ n ] to produce [A'|M], in which the rows of A' again form a quasi- 
orthonormal list of vectors. In this case the solution is 

x p = Gb 

Xh = Py 
P = I n -GA 
G = (A')*M 



where y G J 7 ™. Finally, note that in this variation it is still possible to 
compute P as I n — (A')* A', and to compute x p as (A')* Mb. In other words, 
one doesn't have to actually form G to compute the solution. Finally, while 
this case requires storage for M, it also allows one to compute a solution for 
additional vectors b without having to repeat the ROP step. 

These differences, while trivial for small example problems, may become 
significant when the matrix dimensions are large, and computer storage space 
is limited. Otherwise, choosing one variation over the other is mainly a 
matter of convenience. Pseudocode for these variations is given in Appendix 
B. 



Penrose Conditions 

It is convenient to classify a generalized inverse (G) according to which Pen- 
rose identities it satisifes. In particular, different properties of the solution 
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Gb follow if certain sets of Penrose identities are satisfied. The first Penrose 
identity is AG A = A, which is seen to always be true for our solution. 

AG A = A{A')*MA 

= M- 1 MA(A)*A 
= M~ l A(A)*A 
= M' X X AI A' 
= M- l A' 
= A 

Likewise, the second identity GAG = G is always true. 

GAG = [(A)*M\A[(A)*M\ 
= (A)*A(A)*M 
= (A)*X A ,M 
= (A')*M 
= G 

The third Penrose identity (AG = (AG)*) is seen to be problematic 

AG = A(A')*M 

= M~ 1 MA(A')*M 
= M- l A'(A'fM 
= M- X X A ,M 

If A is of full row rank, then X A > equals I m , and AG becomes I m ; the identity 
is satisfied, allowing the identity to be satisfied. However, if A is not of full 
row rank, then this identity is not true in general. Finally, the fourth Penrose 
identity GA = (GA)* is seen to be true: 

GA = (A)* MA 

= (Ay a 
= ((Ay Ay 

= ((A)* MA)* 
= (GA)* 
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In summary, the generalized inverse obtained by this algorithm is at least 
a {124}-inverse. If it's additionally true that A is of full row rank, then 
it becomes a {1234}-inverse, a.k.a. a Moore-Penrose inverse. Recall that 
generalized inverses which are at least of type {14} (which is the case here) 
yield minimum norm solutions Gb. 

Online Capabilities 

To the author's knowledge, algorithms for solving Ax = b require that all 
data (i.e., A, b) be available at the outset before the linear solver can begin. 
This algorithm is different: it can do the calculation as the data arrives (so 
long as it arrives in a certain manner). This style of computation is referred 
to as an online algorthm. Examples of when it is useful is in cases where 
it takes a large amount of time to compute (or acquire) all the entries in A 
and b. This approach reduces the overall computation time. In addition, 
it will be shown that the updates to x are mutually orthogonal; this means 
that the estimation of \\x\\ is monotonically non-decreasing throughout the 
computation. 

Assume that the rows of [A\b] become available one at a time, and for 
simplicty let the i-th row of [A\b] be the i-th to arrive. Note that when the 
ROP is based on the classical Gram-Schmidt (CGS) 0, E] procedure, the 
1st through i-th rows of A, M and b will no longer change following the i-th 
step of the algorithm. Since those rows are done changing at these points, 
they become available to be used in a computation of x p . The next step is 
to use the column-row expansion [5] on the product (A')*b' to rewrite x p as 



x p = {jvyu = 2^ x f 



(A f yb> = J2 

8=1 

where 

xf = CoL.p')1 K 

and Coh signifies the i-th column. Following the i-th step in the algorithm, 
the i-th term on the right-hand side (ie, x p ) may be computed. Thus the 
solution is accrued just by adding x p + x p + ■ ■ ■ . Furthermore, the up- 
dates to x p are mutually orthogonal, i.e., < x p , , x p >= for i =£ k. This 
approach based on the first variation will be the basis of the following online 
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computation. An illustration of this technique on a numerical example is in 
Appendix D. 

The same basic approach can also be taken for the second variation of 
the algorithm, except now the column-row expansion is used to rewrite the 
generalized inverse as 

m 

G = (A')*M = J2 G(i) 

where 

G» =CoL[(A')*] Rowi[M] 

Note that G^ 1 ' may be computed following the i-th step of the ROP. Following 
the computation of G, the particular solution is easily computed from it. 

Final Remarks 

What is immediately noteworthy about the algorithm [13] is that it doesn't 
use an elimination or partitioning strategy; in particular, there was no solv- 
ing of a triangular system of equations. The algorithm is similar to those 
based on matrix decompositions |10j . in which A is written as a product of 
matrices, and then reinserted into Ax = b. In the algorithm presented here, 
the factorization was done implicitly, by operating on A while it was in the 
context of the equation Ax = b. Borrowing a term from metallurgy, this type 
of factorization might be called an "in situ factorization" . Also, keep in mind 
that it is required that the equations represented by Ax = b be consistent. 
(Although, if they are inconsistent, a simple modification to the pseudocode 
makes it easy to discover that.) 

The new algorithm is perhaps most similar to a version of the QR algo- 
rithm, in which the GS procedure operates on the rows of A. However, in 
that approach, one still has to solve a triangular system of equations. This 
destroys the possibility of easily computing a generalized inverse or a null 
space projection operator, as well as formulating the solution in an online 
mode. 

The solution found by the new method is always of minimum norm. When 
A additionally has full row rank, the solution is a least-square solution. These 
properties follow from the generalized inverse, which is type {124} generally, 
and type {1234} when A has full row rank. (Recall that if a generalized 
inverse G is at least type {14} the solution Gb has minimum norm, and if it 
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is at least type {13} the solution Gb is a least squares solution.) Finally, it's 
pointed out that it isn't necessary to explicitly form G; the first variation 
side-steps that computation. 

The orthonormalization procedure named "ROP" can be thought of as 
just the GS procedure, except that the resulting zero-norm row vectors are 
retained in the end result. Zero vectors are allowed to persist in A only 
because it's easier to leave them there. They could also be removed; in that 
case A and b would have to be re-sized. An extension of the ROP [121 02] 
takes into accout numerical precision and declares the norm of a vector to be 
zero if it is less than a small number e; the size of e is related to the machine 
precision used to implement the algorithm [7J. 

Although the new algorithm was cast to solve Ax = b, it can easily solve 
[12] its matrix generalization: AX = B, where X is n-by-p and B is m-by-p 
(and p > 1). In that case the ROP proceeds as before, and the particular 
and homogeneous parts of the solution X = X p + Xh are 

X„ = GB 
X h = PY 

where G and P are the same as before, and Y G J rnx P, This also admits an 
online formulation. 
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APPENDIX A: pseudocode 

The following pseudocode, written in the style of the C programming lan- 
guage, illustrates the workings of the algorithm. It shows that while it was 
previously expedient to emphasize the role of the matrices M s , it's not nec- 
essary to explicitly form them. Separate pseudocode is presented for each of 
the variations, to aid exposition. 

Two accomodations are made for machine precision, should this be im- 
plemented on a computer. The first is the variable 'eps', which might be 
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set to some multiple of the machine precision (cf. "e-rank" [7]). (For the 
examples herein, 'eps' is zero.) Also, it might be expedient to take further 
action on a row that has a norm less than 'eps'; this would be done in the 
code where the comment "zero- norm option" appears. (However, this option 
is not used in the examples herein.) Finally, the notation RowjP] indicates 
the i-th row of a matrix P. The bars | | indicate a Euclidean norm, and 
the angle brackets <, > indicate an inner product. Since exact arithmetic is 
assumed for the examples and the further discussion in this paper, take 'eps' 
to be zero. 

In the first variation the algorithm begins with the input data A and 
b. Row operations are done on the augmented matrix [A\b], transforming it 
into [A' | b'}. The rows of A' subsequently form a quasi-orthonormal set. This 
version of the ROP is based on the modified Gram-Schmidt (MGS) [6], [TT] 
procedure. 

//first variation 
for (i = 1 to m){ 

mag = | | Row_i [A] I I 
if( mag > eps){ 
//normalization 
Row_i [A] = Row_i [A] / mag 
b_i = b_i / mag 

//orthogonal izat ion 
for (k = i+1 to m){ 

prod = < Row_k[A], Row_i [A] > 
Row_k[A] = Row_k[A] - Row_i [A] * prod 
b_k = b_k - b_i * prod 
} 
}-else-[ 

//implement a "zero-norm option" 
> 
> 

Following this the various solution features (i.e., x p , G, P, ...) are computed. 
In the second variation, the ROP begins with the input of the coeffi- 
cient matrix A and an m-by-m matrix M which is initialized as an identity 
matrix. In the pseudocode, the entries for A and M = I m will be written 
over; the result at the end will be identified as A' and M — (■ ■ ■ M 2 Mi), 
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respectively. The ROP proceeds by doing row operations on the augmented 
matrix [A|7 m ], transforming it into [A'|M], such that the rows of A' form a 
quasi-orthonormal set. This version of the ROP is also based on the MGS 
procedure. 

//second variation 
for (i = 1 to m){ 

mag = | | Row_i [A] I I 
if( mag > eps){ 
//normalization 
Row_i [A] = Row_i [A] / mag 
Row_i [M] = Row_i [M] / mag 

//orthogonal izat ion 
for (k = i+1 to m){ 

prod = < Row_k[A], Row_i [A] > 
Row_k[A] = Row_k[A] - Row_i [A] * prod 
Row_k[M] = Row_k[M] - Row_i [M] * prod 
} 
}else{ 

//implement a "zero-norm option" 
} 
} 

Following this the various solution features are computed. 



APPENDIX B: example of 1st variation 

In this section the pseudocode for the first variation is used to compute the 
solution. The input data are 



A 



-3z 







I 1 


1% 1 


-1 


b = 


2z 


M 2 -Si 


-2 




\l + 4i 



Note that A has rank 2. 
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Step 1 

The % — 1 case in the loop in the pseudocode involves the normalization of 
the first row. The associated row operation is 

1 
Rowi 4— (-)Rowi 

Following this operation, the intermediate values for A' and b' are 



[A\b] 



-i 





i 

3 


2% 1 


-1 


2i 


U 2 — 3i 


-2 


l+4i_ 



Step 2 

This step is the orthogonalization between the first and second rows of A; 
it occurs when i = 1 and k = 2 in the pseudocode. The associated row 
operation is 

R0W2 <— Row 2 — (i)Rowi 

This leads to the following intermediate values for A' and b' 

-i I 



[A\b] 



-i 

2i -1 

H 2 — 3z -2 



ft 

3* 

1 + 42 



Step 3 

This orthogonalization step is between the first and third rows, and occurs 
when i=l and k=3. The associated row operation is 

R0W3 <— R0W3 — (3 + 2i)Rowi 

At this point the intermediate values for A' and b' are 



[A\b] 






—i 





1 " 
3 


2i 





-1 


3 


U 





-2 


3 . 



2 TEE ALGORITHM 



14 



Step 4 

The normalization step for i=2 is for the second row. The associated row 
operation is 

Row 2 ^— (^p)Row 2 
V5 

Following this, A' and b' take on the intermediate values 



[A\b] 






—i 





i 

3 


2 A 





i 


V5„ 


7E l 


VE 


~ l 


4i 





-2 


10 i 
3 ' 



Step 5 

For i=2 and k=3, the third row is orthogonalized with respect to the second. 
The row operation associated with this is 

R0W3 ^— R0W3 — (2v5)Row2 

The intermediate values for A 1 and b' are now 



[A\b] 



-1 







v^ VE I 3 l 

I 



Noteworthy is that the third row of the augmented matrix has become zero. 
Had the third entry in b' been nonzero the equations would be inconsistent. 



Final steps 

According to the pseudocode, the final step should be the normalization of 
the third row (for i=3). However, since the third row has zero norm, this 
step is skipped (according to the if-statement in the pseudocode, since 'eps' 
is zero). It follows that what were identified as intermediate values for A' 
and V in the fifth step are in fact final values. The particular solution may 
now be formed as 



(A')*V 



"0 


-* 


0" 


i 











1 

y/5 
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Also, the null space projection operator is 



P 



(A')*A' 



1 

5 



1 


-2 








1i 


4 



2i 



As noted earlier, the homogeneous solution is formed as Xh = Py, for arbi- 
trary y. Setting y = (?/i, ?/2,2/3) T , where each entry is arbitrary, the result 
is 

i ('\ ( l 

x h = -(2/1 -2iy 3 ) = a 
b \2i \2i. 



where a is an arbitrary element in C. Noteworthy is that the nullity of A 
is one, which corresponds to the above parametrization of Xh requiring only 
one vector. Also note that it was not necessary to explicitly construct the 
M s that were used in the derivation of the algorithm. 



APPENDIX C: example of 2nd variation 

This is the same example as before, the same steps are done, except now M 
is changed instead of b. The row operations Have the same steps and the 
same row operation as in the example for the first variation. Hence, all that 
will be shown are the intermediate vales for the augmented matrix [A|M]. 
Step 1 

0" 
[A\M] = 2% 1 -1 ; 6 1 

1 



Step 2 



Step 3 



\A\M\ 



[A\M] 



"0 -i 


i 

3 


2% 1 -1 


1 o 


4z 2-3z -2 


o 


-i 


1 

3 


2% -1 


-¥ 


4i 2 — 3i -2 





-i | 


1 

3 


2% -1 


-¥ 


a o -2 j - 


- 1 - ¥ 





1 
1 





1 
1 
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Step 4 



[A\M}= jt 



Step 5 



[A\M] 



' 


—i 


o 1 


1 

3 


-2-i 





1 


>/5 


VI I 


15 


_4i 





-2 | 


-1 - 


" 


— z 





1 
3 


*< 





1 


V5 


^5 


15 











-1 




VE 

5 





VE 



The final A' and M matrices are those found above, in step 5. The generalized 
inverse G is found to be 



G = (A')*M 



1 
15 



-2 


— 6i 


yi 





i 


-3 



Besides permitting a computation of the particular solution via x p = Gb, 
it also allows an easy computation of the null space projection operator 
(P = I n — GA), which is the same as before. In short, the second variation 
can compute everything that the first variation did; the difference is that it 
can also compute G. Finally, note that it was not necessary to explicitly 
construct the M s that were used in the derivation of the algorithm. 



Additional remarks 

To complete the example, the matrices M s (s = 1,2,3,4,5) are shown for 
the row operations given in the above examples. (M s corresponds to the row 
operation in the s-th step.) While it was not necessary to actually construct 
and use these matrices, they may provide an aid to understanding for the 
reader. 



Mi 



1 


0" 






1 



Mo 



" 1 


0" 


-i 1 








1 



M, 



1 



-(3 + 2i) 




1 




Ma 



"1 





0" 





V5 
5 





_0 





1_ 



M R 



It is left as an exercise to verify that M5 
puted in the second variation above. 



1 
1 
-2y/E 1 

Mi equals the M that was com- 
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APPENDIX D: online case 

In the following example, the same data is used to illustrate the online case 
for the first variation. In the expressions below, the " input data" for the rows 
of A and b are from newly acquired data. (Double-hyphens in a matrix mean 
that no data has been entered there yet.) Also, the "intermediate results" 
are how A, b, and x p appear following the i-th update. After i — 3, the 
updates are complete, and the last A and b may be identified with A' and b', 
respectively. Also, note how the norm of x p is non- decreasing with respect 
to i. 

i = 1 

input data: 

Rowi(^) = (0,-3i,0) 
Rowi(6) = (1) 



intermediate results: 


A 



—% 



x 



(i) 



input data: 



Row 2 (A) = (2i, 1, 
Row 2 (6) = (2i) 



intermediate results: 



A 



-i ' 
*fii -^ 

5 5 



' I 1 
3 

3 ' 



X 



(2) 



2 ' 

:? 



-¥. 
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input data: 



Row 3 (A) = (4z,2-3i, 
Row 3 (6) = (l + 4i) 



-2) 



intermediate results: 



.4 






— « 


" 




- i - 

3 




"0" 


*#* 





_v / 5 


, b = 


^ 


x (3) _ 





5 




5 




3 


' P 















_ _ 








The particular solution follows from adding all the updates, giving 



•-'■ ■■■■■/•; 



.r 



(1) , (2) , (3) _ i 

T x p T x p — 



which is the same as found earlier. As a final point, note that it is trivial 
to repeat the (CGS) orthonormalization step for each i during this online 
computation. Doing so would increase the accuracy of the solution [12]. 
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